Spontaneous imbibition in disordered porous solids: a theoretical study of helium in 

silica aerogels 



F. Leoni 

GIT- SPEC, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France 

E. Kierlik, M. L. Rosinberg, and G. Tarjus 
Laboratoire de Physique Theorique de la Matiere Condensee, 
Umversite Pierre et Marie Curie, CNRS-UMR 7600 
4 place Jussieu, 75252 Paris Cedex 05, France 

We present a theoretical study of spontaneous imbibition of liquid *He in silica aerogels focusing on 
the effect of porosity on the fluid dynamical behavior. We adopt a coarse-grained three-dimensional 
lattice-gas description like in previous studies of gas adsorption and capillary condensation, and 
use a dynamical mean-field theory, assuming that capillary disorder predominates over permeability 
disorder as in recent phase-field models of spontaneous imbibition. Our results reveal a remarkable 
connection between imbibition and adsorption as also suggested by recent experiments. The im- 
bibition front is always preceded by a precursor film and the classical Lucas- Washburn ^/i scaling 
law is generally recovered, although some deviations may exist at large porosity. Moreover, the 
interface roughening is modified by wetting and confinement effects. Our results suggest that the 
interpretation of the recent experiments should be revised. 

PACS numbers: 47.61.-k,68.15.-fe,47.55.nb,47.56.-|-r 



I. INTRODUCTION 

This paper is devoted to the theoretical description 
of spontaneous imbibition of helium in silica aerogels. 
It is motivated by recent experiments showing a strik- 
ing influence of porosity and temperature on the fluid 
dynamical behavior^. In particular, an intriguing two- 
step process has been observed in very light aerogels at 
low temperature, with two imbibition fronts rising succes- 
sively inside the solid, both obeying the classical Lucas- 
Washburn (LW) linear relationship between penetration 
depth and square-root of tim^. Although these experi- 
ments have not yet been reproduced and cannot be eas- 
ily interpreted at the microscopic scale, they suggest the 
presence of precursor wetting films moving ahead of the 
primary front^, as often observed in other systems, such 
as oil inside glass beads^, water in paper fibers- , or dur- 
ing forced imbibition with immiscible fluid^. More gen- 
erally, these experiments illustrate the subtle interplay 
between disorder, confinement, and wetting phenomena 
that takes place in aerogels and that is also at the origin 
of the behavior observed in adsorption and desorption 
experiment^^t^. Silica aerogels are interesting materials 
for such studies because their porosity can be varied in a 
significant range and one can thus perform a systematic 
study of the influence of the microstructure on the fluid 
behavior. On the other hand, the heterogeneity of the 
gel at many scales and the absence of a well-defined pore 
size make the dynamics more com plex th an in other dis- 
ordered solids, such as Vycor glass'^^'l3_ These features 
also greatly complicate the theoretical description. 

Since spontaneous imbibition plays a crucial role in 
many problems of general interest (from oil recovery to 
agriculture) and is also used as a practical method to 



characterize porous structures (by modeling the solid as 
a bundle of uniform capillaries), it has received con- 
siderable attention over the years^^. In spite of the 
diversity of actual porous media, the ^/i scaling law, 
which results macroscopically from the balance between 
a constant driving capillary force and an increasing vis- 
cous drag, is ubiquitously observed. In fact, it appears 
that the LW relationship holds down to th e nanoscopic 
scale, as shown by recent experimentaP^EH and theoret- 
ical studieJiSHIl! This of course is an important indi- 
cation for the development of nanofluidic devices. The 
dynamics of imbibition fronts is also a subject of much 
interest within the statistical physics community as it 
is a good example of the propagation of an interface 
in a disordered medium (see Ref. [T21 for a re view and 
references therein). Recent experimentaP^'^ni and the- 
oretical invest igat ion^^^'^iiaai j^ave especially focused on 
the roughening of the moving front that results from the 
random spatial distribution of capillary forces. It has 
been shown that the scaling properties of the interface 
are strongly affected by the non-locality induced by fluid 
conservation. Because a microscopic treatment of liquid 
flow in a random medium is a formidable problem, these 
theoretical approaches generally use a coarse-grained de- 
scription in which the spatial conflgurations of the fluid 
are represented by a phase-field modeP^^^ and the disor- 
dered structure of the solid is represented by quenched, 
spatially uncorrelated random fields. While this descrip- 
tion is appropriate to investigate the large-scale prop- 
erties of the interface, it does not account for the cor- 
relations in the solid microstructure (for instance along 
the silica strands in aerogels) nor for the wetting layers 
that may be present on the pore walls, depending on the 
solid wettability. In order to interpret the imbibition ex- 
periments in aeroge one therefore needs a somewhat 
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more involved (though still coarse-grained) description in 
which the solid microstructure is taken into account, and 
the temperature, the porosity, and the solid-fluid inter- 
action can be varied independently. 

In recent years, such a description, based on a lattice- 
gas model and a mean-field statistical treatment, has 
been used to investigate the hysteretic capillary conden- 
sation of gas in disordered porous solids'^"^^^'^. In partic- 
ular, the experimental results in aerogel^^l^ have been 
faithfully reproduced, and the remarkable dependence of 
the adsorption isotherms on temperature and porosity 
has been elucidatecpSHSSI Jn the present paper, we ap- 
ply the same model to spontaneous imbibition, using a 
dynamical version of lattice mean-field theorjEI! that has 
already proven useful to mod el con densation and evap- 
oration processes in nanoporeJ^SEII, Like in the phase- 
field models of imbibitiorPSsSI^ ^ffQ consider the regime 
where capillary disorder predominates over permeabil- 
ity disorder, which corresponds to the case of a slowly 
advancing interface. Only diffusion-like mechanisms are 
thus included and the effect of the hydrodynamic modes 
are incorporated through effective parameters. However, 
the key ingredient at the origin of the constant slowing 
down of the imbibition front predicted by the LW law, 
namely liquid conservation, is explicitly treated. In a pre- 
vious worlP^, we have applied this theory to spontaneous 
imbibition in a single slit-pore, emphasizing the role of 
precursor wetting films on the propagation of the main 
meniscus. This was a first necessary step before address- 
ing the much more complex case of disordered solids. 

The paper is organized as follows. In the next sec- 
tion, we review the lattice gas model and the dynamical 
mean-field theory. The numerical results are presented 
and discussed in Sec. III. Finally, in Sec. IV, we summa- 
rize our findings and conclude. Some details about the 
definition of a local, disorder-dependent permeability are 
given in the Appendix. 



II. MODEL AND THEORY 

As described in detail in Refs. we model the gel- 

fluid system by a coarse-grained nearest-neighbor lattice 
gas with a configurational energy function given by 

(1) 

where = 0, 1 is the fluid occupancy variable at site i 
(i = 1, N) and ?7i = 1, is a quenched random variable 
that describes the solid microstructure {rji = if site i 
is occupied by a gel particle); Wff and Wsf denote the 
fluid-fluid and solid-fiuid attractive interactions, respec- 
tively, and the sums run over all distinct pairs of nearest- 
neighbor (n.n.) sites. The ratio a = Wsf/wff controls 
the wettability of the solid surface and (p = (1/A^) J2i Vi 
defines the solid porosity. 

Gel configurations (i.e., the set of ?7i's) are gener- 



ated by a standard on-lattice diffusion-limited cluster- 
cluster aggregation (DLCA) algorithni'^'^ which has been 
shown to faithfully represent the structure of the base- 
catalyzed silica aerogels which are considered here. A 
body-centered-cubic (bcc) lattice is used in order to avoid 
some lattice artifacts at low temperature^^. From the 
computation of the gel correlation length ^g, the lattice 
constant a (hereafter taken as the unit length) was esti- 
mated in Ref. HH]to be of the order of a few nanometers, 
which is consistent with the coarse-grained picture of a 
gel site representing a Si02 particle with a diameter of 
about 30A. The lattice has lateral size = Ly = L in 
the X and y directions and height Lz in the z direction 
(i.e., N = 2LzL^j^. Periodic boundary conditions are 
taken in the x and y directions and the liquid reservoir 
is located at the bottom of the simulation box, as illus- 
trated in Fig. 1. 




and gel particles are coded in blue and black, respectively, 
and only the fluid sites with a density pi larger than 0.5 are 
shown. The liquid reservoir in located at the bottom of the 
simulation box. 



The theoretical approach that we use to model the liq- 
uid dynamics may be viewed as a mean-field version of 
the Kawasaki spin exchange dynamics and was originally 
applied to study phase separation and surface enrichment 
in solid binary alloy^^ (see Ref. dHl for a comprehen- 
sive review, Ref. |36| for an application to diffusion in 
membranes and bulk fluids and Refs. [301 and [31] for a 
recent application to fluids confined in nanopores). It 
is referred to as "mean-field kinetic theory" or "dynamic 
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mean- field theory " (DMFT) . The starting point is the ex- 
act master equation describing the temporal evolution of 
Pi{t) = {TiTjiit, the ensemble average fluid density at site 
i. Assuming that transport is due to a hopping process 
between nearest-neighbor sites and replacing the occu- 
pancy variables by their ensemble average (z. e. neglecting 
dynamical correlations) , one can then derive a continuity 
equation describing the local conservation of mass 



dt 



j/i 



(2) 



where the sum runs over all n.n. of site i and the flux 
Jij{t) is given by 



Jij{t) = w,jPi{rij - Pj)- w-jiPjifj, - Pj) 



(3) 



(note the presence of the random variables rji and rjj that 
describe the configuration of the gel particles on the lat- 
tice). The Wij's are the Metropolis transition probabil- 
ities associated to the Kawasaki dynamics in the mean 
field approximation, 



where 



and 



0, 

Ej — Ei, 



Ej < Ei 
Ei > Ei 



j7» 



(4) 



(5) 



(6) 



as obtained from Eq. (1); wq is an elementary jump rate 
that sets the time scale and is considered here as an effec- 
tive parameter incorporating some information about the 
hydrodynamic modes that are not explicitly treated. It 
is indeed clear that this approximate treatment of the ac- 
tual dynamics is essentially diffusive, in the same way as 
the recent mesoscopic phase-field models of spontaneous 
and forced-flow imbibitior t^^l^^T^ . The analogy with a 
continuum phase-field description becomes more appar- 
ent if the continuity equation is recast into the discrete 
version of a generalized Cahn-Hill iard equatiorP^, or a 
dynamic density functional theorji2S129|^ introducing an 
effective mobility coefficient related to the local site den- 
sities via the Metropolis transition probabilities'^" (see in 
particular the appendix of Ref . |32] where the continuum 
limit of the above equations is derived in the much sim- 
pler case of a single slit pore). 

The connection of Eq. ^ to the phenomenological 
Darcy's equatiorP^l that relates the fiux of liquid in a 
porous medium to the macroscopic pressure gradient. 



J = --Vp 



(7) 



is established by taking proportional to the ratio k/t] 
of the solid permeability to the fluid viscosity (as in the 



phase-field models^. The permeability k represents the 
(volume-averaged) frictional effects exerted by the solid 
microstructure on the imbibing fluid and is essentially de- 
pendent on the size of the pores through which the liquid 
flows (for a single capillary of radius R, k (x if one as- 
sumes a Hagen-Poiseuille laminar flow) . For porous solids 
with a reasonable homogeneous cylindrical pore structure 
such as Vycor, k can thus be related to phenomenologi- 
cal parameters characterizing the solid morphology such 
as the p orosit y (j), the average pore radius rp, or the tor- 
tuosity TEnmi. This relationship is much more elusive in 
aerogels since the structure is inhomogeneous at many 
scales, as already emphasized. However, in order to com- 
pare results obtained in aerogels with different porosi- 
ties, it is crucial to take into account some dependence 
of Wo on the solid properties. In a first approximation, 
one may simply assume that k, and consequently wq, 
only depends on (f), neglecting local variations due to 
the disorder of the gel structure. Several theoretical ex- 
pressions relating k to (p are available in the literatur^^. 
In particular, Brinkman's effective-medium theorjPl' ap- 
pears to be a good representation of the flow in a di- 
lute, disordered porous matrix down to the nanoscale, as 
shown by very recent lattice-Boltzmann and molecular 
simulations^. For a dilute collection of non-overlapping 
spheres of radius R and volume density n (so that the vol- 
ume fraction is l — (j) = AnR'^n/S), Brinkman's expression 
for k(0) is 



k(0) = Kq 



1 



(8) 



where kq — l/(67ri?n) — 2i?^/[9(l — (j))] is the perme- 
ability in the high-dilution limitp2l. Brinkman's theory 
somewhat underestimates the permeability of a dilute 
gel network'*^ but the correction is almost independent 
of (f) in the regime considered in the present study 
{(j> > 87%). The transposition to the case of a discrete 
b.c.c. lattice suggests to take R/a — ■\/3/4 so that 
two spheres representing silica particles and placed 
at n.n. sites are tangent (this choice, however, is not 
essential in our treatment since the value of R is the 
same for all aerogels). This yields the values indicated 
in Table 1. Overall, assuming that wq has the same 
dependence on (p as the permeability simply amounts to 
replacing the imbibition time t in Eq. ([2| by a rescaled, 
porosity-dependent, time t(0) — t/K{(j)). 








0.87 


0.113 


0.90 


0.174 


0.92 


0.244 


0.95 


0.473 



TABLE I: Permeability K((t>) computed from Eq. ^ 
R/a = 73/4. 



with 
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Of course, this is only an average description and one 
would like to also take into account the influence of lo- 
cal inhomogeneities in the microstructure. This can be 
done for instance by defining a local porosity obtained 
by averaging the number of gel particles in a spherical 
region around each fluid site. Such an approximate treat- 
ment, which introduces some disorder into the permeabil- 
ity (hence into wq), is described in the Appendix. It turns 
out, however, that this correction is almost negligible. In- 
deed, as emphasized in Refs. I12I21I in the framework of 
phase-field models, the effect of disorder on the mobil- 
ity is essentially proportional to the interfacial velocity, 
which is never large in the case of spontaneous imbibition 
and moreover continuously slows down (as t~^^'^ accord- 
ing to the macroscopic LW law). Therefore, disorder- 
induced fluctuations of the capillary forces should prevail 
over fluctuations of the permeability, and the former are 
properly taken into account in our treatment via the de- 
pendence of local energy Ei on the quenched variables rji 
representing the solid microstructure. 

Finally, it is important to observe that the local den- 
sities {pi} minimizing the static (mean- field) Helmholtz 
free energy J'^[{pi}] in the canonical ensemble 

/?-^[{Pi}] Y^iPi ln(pO + iVi - Pi) ln(??j - pO] 

i 

- Pwff ^ p^Pj - PWsf ^ [p^{l - Tlj) + Pjil - Vi)] 
<ij> <ij> 

(9) 

(with /3 — l/ksT) or the static (mean- field) grand- 
potential fl = T — p J2i Pi in t-iis grand canonical ensem- 
ble, are steady-state solutions of the DMFT equation^^S] 
Indeed, one can readily check that the pi's satisfying the 
set of coupled equations dT/dpi = p, i.e. 



Vi - Pi 



ejip[l3{p ~ Ei)] 



(10) 



are solutions of Eq. ^ when dpi{t)/dt ~ 0. This feature 
will allow us to unveil a remarkable connection between 
the microscopic states visited dynamically during the im- 
bibition process and those visited during quasi-static ad- 
sorption experiments performed in the same systems. In 
a disordered material, Eq. ( 10 1 has many solutions in 



general, typically of the order of exp{N), and these local 
minima^* correspond to metastable states whose pres- 
ence underlies the out-of-equilibrium behavior observed 
in adsorption/desorption cxperimcntd^SttSSj^ This feature 
also plays a crucial role in imbibition, as the interface 
may be pinned in one of these states, resulting in an in- 
termittent, avalanche-like motiorfi^l. 

In practice, we have integrated Eq. ^ via Euler's 
method [i.e., pi{t + At) = Pt{t)-AtJ2j/i Jij{t)] by using 
the dimensionless time step u;o((/>)At = 0.1. This value 
was small enough to ensure the stability of the solution in 
all cases. However, it corresponds to a very small varia- 
tion of the local site densities. This was already found in 



our previous study of imbibition in a single slit-por^^Sl^ 
but the present fully three-dimensional computation is 
considerably more time-consuming. As a result, the pen- 
etration depth of the main interface remains small com- 
pared to the overall height of the simulation box, even 
after several months of CPU time. This is a serious tech- 
nical problem which, unfortunately, seems to be inher- 
ent to this type of approach. Despite a great deal of 
effort, we have not found an alternative method. For 
plotting the numerical results, we use a larger time scale, 
Iq — W^woAt, which is hereafter taken as the effective 
time unit. 



III. RESULTS AND DISCUSSION 

We have studied aerogels with porosity (j) = 87%, 90%, 
92%, and 95% but for the sake of brevity we mainly focus 
on the 87% and 95% solids in the following. The results 
correspond to the case of a primary imbibition, when 
the liquid in the reservoir is initially in contact with the 
porous solid filled by (bulk) vapor at p = Psa^ (this is 
not an equilibrium condition due to the presence of the 
porous solid) . This is realized by fixing the density of the 
fluid at the saturated liquid density p; (T) in the first layer 
at the bottom of the box (z — 0) and at the saturated 
vapor density Pg{T) in the rest of the gel. Then, at t = 0, 
we let the system evolve according to Eq. ^ and the 
liquid enters the gel. 

To make contact with previous theoretical studies of 
"'He adsorption^-^, the interaction ratio a = Wgf/wff 
is set equal to 2 in order to reproduce approximately 
the height of the sorption hysteresis loop measured in 
a 87% aerogel at T = 2.42K. However, to prevent the 
liquid from escaping, a drying condition is imposed in 
the top layer by setting a = for z = L^. All cal- 
culations are performed at the same reduced tempera- 
ture T* = kBT/wff = 0.9 {i.e. T/Tc = 0.45) so that 
Pg « 0.013 and pi « 0.987. 

As emphasized in Ref. |32l the possible presence of 
precursor films that advance faster than the main inter- 
face and that quickly reach the other end of the porous 
solid makes the system size Lz in the direction of the 
flow a critical parameter in this type of studj^Sl. This 
is illustrated in Fig. [2] where we plot the time evolution 
of the fluid density profile p(z, t) obtained by averaging 
the density in the x — y plane and we compare the re- 
sults for two samples of a 95% porosity aerogel with sizes 
Lz = 2L = 200 and = 4L = 400. The characteris- 
tic features of the profile will be discussed in more detail 
below but one can see at once that because of the bound- 
ary condition imposed &t z = Lz the fluid accumulates 
at the top end of the solid, which induces an artificial 
thickening of the density profile (this is not due to a "re- 
sistance" of the gas occupying the gel during the liquid 
invasion^^ since the gas can always condense). However, 
since it takes some time for the precursor film to reach 
the other end of the solid (see below), these finite-size ef- 
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FIG. 2; (Color on line) Profiles of the average fluid density 
p{z, t) at different imbibition times {t = 30, 60, 90, 120 from 
bottom to top) in two 95% aerogels with sizes = 2L = 200 
(red lines) and = — 400 (blue lines). The latter sample 
has been duplicated from the formei^^. The length and time 
units are o and to, respectively. 
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fects can be neglected so long as the imbibition time t is 
shorter than some time tmax{z) (the larger z, the smaller 
irnax)- For instance, in Fig. [2j the profiles obtained with 
— 200 are insensitive to the boundary effects if the 
study is restricted to z < 100 and t < 120. If one is only 
interested in the small-z region, typically z < 20, larger 
times may be considered, up to t « 200. The results 
that we discuss hereafter were obtained with a system of 
height L,^=2L = 200 (whence iV = 4 x 10^). 

To begin with, we show in Fig. 3 the time evolution 
of the profiles p{z, t) in the different aerogels, after aver- 
aging the numerical data over a small number of gel re- 
alizations to somewhat smooth out the fluctuations over 
disorder (it was not possible to consider more realiza- 
tions because of the huge computational time required). 
Note that the time unit is to (just proportional to the 
number of time steps in the integration of Eq. ^ by 
Euler's method) and is therefore not renormalized by 
the porosity-dependent factor k(0) given by Eq. ([s]). 
We are indeed only interested here in observing quali- 
tatively the influence of porosity on the density profiles 
and not in comparing the imbibition rates (as will be 
shown below, the imbibition front actually moves faster 
as the porosity increases when time is properly rescaled) . 
The most salient feature in all cases is that the density is 
much larger than pg far ahead of the main interface, in- 
dicating the presence of a precursor "wetting" film that 
thickens with time. Whereas the imbibition front can- 
not be clearly discriminated from the precursor film for 
(f) = 87%, there are clearly two distinct regions in the 
profile for = 95%. In this latter case, the arrival of 
the imbibition front at a given distance z from the liquid 
reservoir is signaled by a steep increase of the density 
from about p « 0.30 to p « p,(l - 0) = 0.94. Note that 
the precursor has already reached the middle of the sam- 



FIG. 3: Density profiles p{z,t) in aerogels with porosity 87%, 
90%, 92% and 95% at times t = 1, 2, 30, 60, 90, 120.. .210. The 
data are averaged over 3 gel realizations for = 87% and 5 
realizations for the other porosities. 



pie at t w 30 while the front is still located at z w 10 (see 
Fig. 3). 

In order to better visualize the underlying microscopic 
mechanisms, we show in Fig. [4] some typical snapshots 
of the liquid configurations in two samples of 87% and 
95% porosity. These snapshots are taken at successive 
times in a vertical plane located in the middle of the 
solid (x = L/2 = 50). 

One can see that the fluid density distribution is es- 
sentially bimodal at this temperature {i.e., pi ss Pg or 
Pi ~ Pi) with a negligible fraction of sites with an in- 
termediate density. This feature was also observed in 
the study of helium adsorption in aerogel^^H i^ fact, 
the configurations of the imbibed liquid closely resem- 
ble those observed during adsorption as one slowly varies 
the chemical potential (or the pressure) in the external 
reservoiil^sHlll This is even more visible in Figs. [5]and[6] 
that display horizontal cross sections of the system at two 
distances z from the liquid reservoir (which is somewhat 
equivalent to examine the same cross section at succes- 
sive times: in this sense, Fig. [6] is a zoom on the first 
stages of Fig. [5]). Like in adsorption, the first stage of 
the imbibition process is a coating of the silica strands 
by a liquid film. (In fact, as soon as the system is let 
to evolve, there is a quasi-instantaneous redistribution of 
the gas molecules due to the attraction exerted by the 
solid, as can be seen for instance in Fig. [6]for i = 2; this 
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FIG. 4: (Color on line) Snapshots of the liquid configurations 
in 87% (top) and 95% (bottom) aerogel samples taken in the 
vertical plane 2; = L/2 = 50 (only the region 2: < L is shown) 
at different times. The aerogel particles are shown in black 
and the fluid density is represented with the color scale in- 
dicated in the figure (from pi ~ in grey to p; = f in dark 
blue). The red lines indicate the position of the horizontal 
planes corresponding to the cross sections shown in Figs. [5] 
andlD 



transient process determines the true initial configuration 
of the fluid inside the porous solid and has noth ing to do 
with the arrival of the liquid from the reservoip22l.) The 
wetting film then progressively thickens while the voids 
or crevices between the aerogel strands fill with liquid. 
During this second stage, the porosity of the solid plays 
a crucial role since the size and shape of the voids are 
quite different in the two aerogels^^ . This has also an 
important consequence on the morphology and the ad- 
vance of the imbibition front, as can be seen in Fig. [4] 
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FIG. 5: (Color on line) Horizontal cross sections in the plane 
2 = 15 at different times t. The samples are the same as in 

Fig.g 



When looking at the density profiles in the 87% and 
95% aerogels in Fig. 3, one cannot help but notice a 
similarity with the corresponding adsorption isotherms 
computed in Refs. I26I27I Those were obtained from 
Eqs. ( 10 1 by starting from a low chemical potential and 



(Color on line) Same as Fig. |5]for z = 50. 



predicts a genuine jump in the average density at low 
enough temperature). To check if there exists a more 
quantitative relationship between the two processes, it 
is instructive to compute the (dynamical) chemical po- 
tential profile /i(z, t) defined as the average of the local 
chemical potential (f) = kBT\n[pi{t)/{l—pi{t))]+Ei{t) 
over the transverse directions. Some profiles fi{z) at dif- 
ferent imbibition times are shown in Fig. [Tj 




increasing /i in small steps. As also observed experimen- 
tally, these isotherms are smooth in the 87% aerogel but 
become much steeper in the 95% one (the theory actually 



FIG. 7: Chemical potential profile p(z, t) in the 87% and 95% 
aerogels at t = 1, 30, 60, 90, 120 (from left to right). 



It is a priori unclear whether this quantity has a true 
physical meaning but it turns out that the /i^'s at a given 
z fluctuate very little from site to site, as illustrated in 
Fig. |8] by a series of histograms collected at different 
times. One can see in Fig. [7] that fi{z, t) ahead of the 
imbibition front is initially much lower than psat due to 
the attraction exerted by the solid (this corresponds to 
the initial redistribution of the vapor inside the aerogel 
noticed above). The effective chemical potential then 
grows with time as the density of the wetting fllm and 
the gas density in the largest cavities both increase. 

Now, plotting p{z,t) as a function of ^isat — ^J'{z,t), as 
done in Fig. |9j reveals that the density proflles at dif- 
ferent imbibition times collapse on a single master curve 
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FIG. 8: (Color on line) Histograms of the local chemical 
potentials fii in the 87% and 95% aerogels at z = 50 and 
t = 30, 60, 90, 120, 150, 180 (from left to right). The small dis- 
persion around the average fi{z, t) over the transverse direc- 
tions shows that the fluid has reached approximately a local 
(metastable) equilibrium in the plane x — y. 



FIG. 9: (Color on line) Dynamic density proflles p{z, t) in the 
87% (a) and 95% (b) aerogels at different times as a function 
of the chemical potentials fi{z,t) (red points): comparison 
with the corresponding gas adsorption isotherms p(p) (black 
solid curves). Note that the fluid density is here plotted as 
a function of Usat — fJ-iz, t) (resp. ^isat — n) instead of fj,{z, t) 
(resp. /i). 



which is very close (though not fully identical) to the cor- 
responding adsorption isotherm. This surprising finding 
strongly suggests that the same fluid configurations (or 
at least very similar configurations) are encountered dur- 
ing the imbibition and adsorption processes. This can be 
checked by comparing the configurations obtained dy- 
namically at time t to those obtained adiabatically in 
the same plane z via a grand-canonical adsorption proto- 
col with fi = i.j,(z,t). This comparison (not shown here) 
proves unambiguously that the same states are indeed 
visited in the 87% aerogel whereas some differences oc- 
cur in the 95% on^^ [whence the small but noteworthy 
differences between the two curves in the vertical part of 
the isotherm plotted in Fig. [o] (b)]. 

How can one rationalize this remarkable result ? 
Firstly we note that the quasi-uniformity of the /ij 's in the 
direction perpendicular to the fiow**^ implies that the flux 
Jij (t) from a site « to a nearest-neighbor site j in the same 
plane is approximately zero (indeed, inserting the defini- 
tion of /ii in Eqs. (4) and (5) allows one to rewrite Eq. 
^ as J.jit) = w;y/9,(l-pj)[l-exp(-/3(/i,-^,))pSI). in 
other words, the fluid "equilibrates" much more rapidly 
in X — y plane than in the direction of the flow (here, 
"equilibrates" means that some local minimum in the 
energy landscape has been reached). This becomes more 
and more accurate as the imbibition front advances and 
slows down, as illustrated by the decrease of the variance 
in the histograms of Fig. [8j Secondly, as can be seen 
in Fig. [7j fj,{z) varies very slowly with z and remains 
essentially constant over all characteristic lengths in the 
system (in particular, the aerogel correlation length ^g; 
see the discussion below). Therefore, at time t, the fluid 
structure at the distance z (or, say, in a slice "around" 



z) is the one observed in the same slice of a three- 
dimensional metastable state, i.e. a local minimum of 
the grand potential Q for fi = ^{z,t). There are actu- 
ally many metastable states for a given value of /i (vis- 
ited when performing incomplete adsorption/desorption 
cyclef|2212£l) ^ but since n{z,t) inside the aerogel is very 
negative at the beginning of the imbibition process, it 
is not surprising that the states visited successively dur- 
ing the process are (to a good approximation) the same 
extremal states as in an adsorption experiment. As a fur- 
ther evidence of this assertion, we compare in Fig. |10| the 
imbibition profile at time t = 60 with the (static) density 
profile obtained by solving Eqs. ( [lO| with the fictitious 
nonuniform external potential Vextit, z) = figat — IJ-{t, z), 
starting from the initial conditions Pi = 0, Vi. Apart from 
small deviations in the region of the imbibition front, one 
observes a remarkable agreement between the two curves. 

The above reasoning, however, only applies when the 
adsorption isotherm is smooth in the thermodynamic 
limit. This is always the case in the 87% porosity aerogel 
because the available empty space between the silica par- 
ticles is small so that the elementary condensation events 
(the so-called "avalanches") remain localized^'' On 
the other hand, in very light aerogels and at low tem- 
perature {e.g., 4> = 95% and kBT/wff = 0.9), a genuine 
discontinuity is predicted to occur in the isotherm in the 
thermodynamic limit at a certain value of /J^. This jump 
is due to the condensation of the gas inside a large void 
space spanning the whole porous solid (in other words, 
the occurrence of a macroscopic avalanche). The nature 
and properties of this nonequilibrium, disorder-induced 



FIG. 10: (Color on line) Comparison between the imbibition 
profile in a 87% porosity aerogel at f = 60 (solid black line) 
and the density profile obtained from Eqs. ( 10 1 using the fic- 
titious external potential Vextit, z) = Usat — fi(t, z) and initial 
conditions Pi = (red dashed line). 



phase transition are discussed in detail in Ref. HSl The 
presence of a discontinuity in the isotherm also signals 
that there are no available metastable states (within a 
certain range of density) in which the system can locahy 
relax. Therefore, within this range of density, the fluid 
configurations visited during the imbibition process can- 
not be genuine metastable states. This explains the de- 
viations observed in Fig. |9](b) in the vicinity of the jump 
in the isotherm, which also corresponds to the range of 
densities associated to the imbibition front. On the other 
hand, there is good agreement between the two curves for 
0.3. 

In the preceding discussion, for instance in the para- 
metric plot p{z, t) vs. iJ,{z, t) shown in Fig. |9j t dropped 
out of the description. We now consider the time evolu- 
tion of the imbibition profiles. Since imbibition proceeds 
in several stages that can be more or less distinguished 
depending on the gel porosity (advance of the precursor 
film, swelling of the film and filling of the voids between 
the silica strands, advance of the main front), it is in- 
structive to study these stages separately although they 
take place simultaneously in the different regions of the 
solid. 

In Fig. [11] we show the time evolution of the aver- 
age density p(z, t) at different heights above the reser- 
voir before the arrival of the main front. We first es- 
timate the time tpre{z) that it takes for p{z,t) to in- 
crease by 10% from its initial gas-like value, which we 
(rather arbitrarily) choose as the signature of the ar- 
rival of the precursor film. In both aerogels, this time 
varies like z'^ (for instance tpre{z) ~ z^/1000 in the 95% 
aerogel) indicating that the rise of the wetting film is 
purely diffusive. This behavior was already observed in 
the slit pord^ in agreement with molecular and lattice- 
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FIG. 11: (Color on line) Time evolution of the average fluid 
density p{z, t) at different heights z in aerogels of 87% (left) 
and 95% (right) porosity (the results have been averaged re- 
spectively over three and five gel realizations). The dashed 
lines correspond to the logarithmic fit p{z, t) = a+bhi[z /t^^'^) 
with a ^ 0.669 (resp. 0.334) and b « -0.220 (resp. -0.10) 
for the 87% (resp. 95%) aerogel. 



Boltzmann simulationJl^IIil^ and is also found in several 
experiments''^^. The density p(z, t) then slowly increases 
with time as the wetting film along the silica strands 
thickens, the gas density within the largest voids in- 
creases and the liquid invades t he s mall crevices between 
the strandd^. As shown in Fig. 11 the numerical data in 
this intermediate time regime are well fitted by a logarith- 
mic growth law (although a power-law dependence with 
a very small exponent cannot be discarded) . We have no 
definite explanation for this behavior but it may simply 
be due to the use of a nearest-neighbor solid-fluid inter- 
action. It is indeed reminiscent of the logarithmic growth 
of wetting layers that is found in the case of short- range 
{i.e. exponentially decreasing) interface potentia P^I^ . 
More important in the present context is the observation 
that p{z,t) « a -|- 61n(z/i^/^) with a and 6 only depend- 
ing on the porosity. This shows that the growth of the 
invading liquid domain follows a purely diffusive scaling 
law in this time regime, despite the complexity of the gel 
microstructure. 



Eventually, the main imbibition front arrives at the 
height z and the average density increases up to pi(l — </)), 
either smoothly for (j) = 87% or more abruptly for 
(/> = 95%. However, even in the lightest aerogel, the 
front is not sharp. Therefore, to analyze its time evolu- 
tion, we extract from the density profiles in Fig. 3 a series 
of rise level curves corresponding to the different heights 
Zp{t) at which the density p is reached for the first time 
(with 0.5 < p < 0.8 for = 87% and 0.4 < p < 0.8 for 
(/> = 95%). As shown in Fig. 12, these curves are almost 
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linear when plotted against t^/^, although the data are 
somewhat noisy due to the small number of gel realiza- 
tions. Accordingly, a reasonable collapse of the imbibi- 
tion fronts at different times is obtained in Fig. [13] by 
using the scahng variable z/t^/^. 




FIG. 12: Rise level curves in the 87% and 95% aerogels as 
a function of t^^'^ . These curves represent the time evolu- 
tion of the position Zp{t) of the imbibition front defined by 
p(zp{t),t) = p. 



average position of the front which we define as 

hit) = — ^ — r^^^ zp{t)dp (11) 

Praax Pmin Jp„ii„ 

with Pmin — 0.5 and Pmax — 0.8 (these values are not 
critical as long as Pmin is large enough to leave aside 
the contribution of the precursor films). It appears that 
a better fit of the data for (j) = 95% is obtained with 
h{t) ~ t'^-*^ (and indeed a slightly better collapse is also 
obtained in Fig. 13 by using the scaling variable z/t^-^^). 
Note that the imbibition times in Fig. 14 are renormal- 
ized by the factor k((/)) to describe the variation of the 
permeability with porosity. 
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FIG. 13: Average density profiles p{z, t) in the 87% and 95% 
aerogels plotted as a function of the scaling variable z/t^^'^. 
The figure focuses on the imbibition front in the time interval 
m<t< 150. 



From this analysis, it is tempting to conclude that 
the advance of the imbibition front obeys the Lucas- 
Washburn ^/t law for all porosities. However, when look- 
ing more closely, it seems that the advance of the front 
is a bit slower in the 95% porosity aerogel. This is illus- 
trated in Fig. 



14 that shows the time evolution of the 



FIG. 14: (Color on line) Time evolutio n of the position h{t) 
of the imbibition front defined by Eq. (Ill with pmin ~ 0.5 
and Pmax = 0.8. The red dashed lines represent the best fit to 
the power law h{t) ~ t" with a = 0.49 for (j) = 87% (bottom 
curve) and a = 0.46 for <j!> — 95% (top curve). The imbibition 
times are renormalized here by the porosity-dependent factor 
it{4>) given in Table 1. 

Considering the limited amount of data at our disposal 
(due to the very slow progression of the imbibition front 
in the numerical calculations), it is unclear whether this 
sub-difFusive behavior actually persists at longer times. 
However, this may be a genuine signature of the fractal- 
ity of the aerogel. Although 95% aerogels displays fractal 
character on a rather limited length scale (and moreover 
the void space filled by the liquid is not fractal) , it is plau- 
sible that the progression of the imbibition front, which 
corresponds to the filling of the largest cavities in the 
gel, is actually sensible to this characteristic feature of 
the microstructure. However, to settle this delicate is- 
sue, it would be necessary to study larger samples and 
even lighter aerogels, which is computationally too de- 
manding at the moment. As shown in Fig 15, the total 
fluid uptake T{t) — J[p{z,t) — p{z,0)]dz shows a similar 
behavior although the exponent a is a bit larger as the 
contribution of the precursor films (which behave purely 
diffusely in all aerogels) is included. 

The results displayed in Figs. 12 and 13 also imply 
that the effective width W{t) of the imbibition front in- 
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FIG. 15: (Color on line) Time evolution of the total fluid up- 
take r{t). From bottom to top: (j> = 87%, 90%, 92%, and 
95%. The exponents of the best power-law fit (red dashed 
lines) are 0.50, 0.49, 0.49 and 0.47, respectively. The imbibi- 
tion times are renormalized by the porosity-dependent factor 
K((ji) given in Table 1. 



creases as a power of time with an exponent /3 close 
to 1/2 (and no smaller than 0.45 for cf) = 95%). Care 
must be exerted in drawing definite conclusions from 
the present calculations^ but they are clearly in dis- 
agreement with the behavior predicted by phase-field 
calculation^^!] -^yjiich suggest /? « 0.3. In this respect, it 
must be emphasized that the disordered matrix in phase- 
field models is described by quenched, spatially uncorre- 
lated random fields. At this level of coarse-graining, there 
is no length scale associated to the solid microstructure 
(for instance, an average pore size Rp) and the growth 
of interface fluctuations is only limited by L, the lateral 
size of the system. These fluctuations are then correlated 
up to a distance ^ which turns out to be directly related 
to the progression of the interface so that it grows with 
time like t^^^. Since the interface is found to be super- 
rough with a global roughness exponent x ~ 1.25, this 
eventually leads to W{t) ~ t^-"^. It is not clear, however, 
whether this line of arguments remains valid when the 
disorder is strong and even the largest pores in the solid 
are small (in a 87% DLCA aerogel for instance, there are 
essentially no cavities with a size larger than 5a, where 
a is the lattice unit; the gel correlation length ^g, which 
is closely related to the average size of the cavities, is of 
the same oide^^. In such a situation, the actual corre- 
lation length in the confined liquid saturates at a value 
corresponding to the size of the cavitieJ^ and one rather 
expect the "trivial" value f3 = 1/2, as if the solid were 
just a bundle of independent pores. In fact, this is the 
value found in recent experiments performed with wa- 
ter and hydrocarbons in porous Vycor glas^^l^ a result 
that was interpreted as signaling the absence of network 
effects during the imbibition process. Confinement ef- 
fects are indeed strong in Vycor which has a low porosity 
{(p « 30%) and a small average pore radius {Rp w 3 — 5 



nm). It is possible, but remains at the speculative level, 
that there is a critical disorder (in the present case, a crit- 
ical value of the porosity) separating a "strong-disorder" 
regime dominated by confinement, with trivial /3 — 1/2 
and a "weak-disorder" regime where the asymptotic (5 is 
close to 0.3. 

Finally, it remains to compare the time evolution of 
the imbibition process in the different aerogels. This re- 
quires to renormalize the time for taking into account the 
influence of the porosity on the permeability K{(j)), as dis- 
cussed in section 2. The results displayed in Fig. 14 show 
that the imbibition front rises faster in the lightest aero- 
gel, which is indeed what is observed in experiment^. 
The same behavior is found for the total fluid uptake 
r{t) displayed in Fig. 15. As discussed in the Appendix, 
these results are essentially unchanged when taking into 
account local variations of the permeability. 



IV. CONCLUSION 

In this paper, we have presented a theoretical study of 
spontaneous imbibition of helium in silica aerogels using 
a coarse-grained lattice model of the gel microstructure 
and a mean- field treatment of the fluid dynamics. Our 
main approximation, which should be valid in the case of 
a slowly advancing front, is that capillary disorder pre- 
dominates over permeability disorder, as also assumed 
in phase-field treatments of spontaneous imbibition. To 
our knowledge, this work is the first three-dimensional 
study of imbibition in a disordered porous solid. The 
main findings can be summarized as follows: 

1) Irrespective of porosity, the first stage of imbibition 
corresponds to the advance of a liquid film along the silica 
strands and in the small crevices of the microstructure. 
The gas density then increases in the larger voids between 
the strands while the precursor film thickens. The largest 
cavities in the aerogel eventually fill with liquid, which 
may be associated to the passage of the imbibition front. 
It is however difficult to clearly distinguish the front from 
the precursor film except in the case of the lightest gels 
{e.g., a 95% porosity aerogel). 

2) The time evolution of the whole process is in general 
purely diffusive but there are some indications that the 
advance of the front could be slower than y/i in very light 
aerogels, which may be related to the fractal character of 
the solid. When taking into account the dependence of 
the 'effective' permeability on the porosity, one finds that 
the front advances faster in light aerogels, as observed 
experiment ally. 

3) The effective width of the interface also varies ap- 
proximately like v^, which is not the behavior predicted 
by phase-field models with uncorrelated random fields. 
Although our results are too limited to draw definite 
conclusions, this indicates that wetting and confinement 
modify the interface roughening. 

4) Due to the quasi-static behavior of the imbibition 
process, the fiuid configurations observed during imbibi- 
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tion and gas adsorption bear a strong resemblance. As 
a consequence, when parametrically plotting the average 
fluid density against the average chemical potential at a 
fixed height and at various imbibition times, one obtains 
a master curve that almost identifies with an adsorp- 
tion isotherm computed in the grand canonical ensem- 
ble. This correspondance breaks down in high-porosity 
aerogels at low temperature when a jump discontinuity is 
present in the adsorption isotherm, signaling a first-order, 
nonequilibrium, disorder-induced phase transition. 

These findings are in partial agreement with the exper- 
imental observations of Ref. tU with aerogels of porosities 
90.4, 95.8 and 99.5%. There is a clear change in the fluid 
behavior as the porosity increases, which is induced by 
the presence of very large cavities. This can be put in 
parallel with the change of behavior observed in adsorp- 
tion isotherms, as correctly suggested by the authors of 
Ref. [H Moreover, it is plausible that the distinction be- 
tween the precursor film and the main imbibition front is 
enhanced visually by the non-linearity of the eye vision, 
and can be interpreted as the existence of two distinct 
advancing interfaces. However, even in the case of the 
lightest aerogel, we find that the first interface is always 
advancing faster than the main front, contrary to what 
has been observed experimentally. We thus rather believe 
that the two advancing interfaces seen in the experiments 
correspond to either a unique interface viewed from dif- 
ferent directions or to two distinct imbibition fronts gen- 
erated by macroscopic inhomogeneities in the aerogel. 
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the location of the first minimum of the two-point cor- 
relation function) but, to illustrate the effect of ro, we 
report in Fig. [16] the histograms of the averaged solid 
densities 1 — (f>i{ro) for several values of rg : for each site 
j, we compute 4>i{rQ) = J2j Vj / J2j 1 where the sum is 
performed over sites j whose distance from site i is less 
than or equal to tq; then, one builds histograms from all 
the empty sites. For aerogels of porosity (j) = 87% and 
4> = 95%, one has — 4a and — 10a, respectiveljffH. 
As expected, when tq is too small compared to ^g, there 
is a non-zero probability that the coarse-grained density 
is zero, which is not acceptable since the corresponding 
Brinkman's permeability then diverges (see Eq. Al I. On 



the other hand, when tq becomes larger, the distribution 
is more peaked around the mean permeability of the ma- 
terial. Taking rg equal to avoids the spurious effect of 
a local permeability that diverges and captures reason- 
ably well the local morphology (i.e., the correlations) of 
the medium. 
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Appendix A: Using a local permeability 

In this Appendix we propose and test a modification 
of the theoretical treatment that takes into account the 
local environment of the fiuid particles inside the gel in 
order to introduce some disorder in the permeability k. 
The idea is to define a local porosity at each site by tak- 
ing into account the effect of the environment within a 
given distance. More specifically, we introduced a coarse- 
grained porosity field 4'iij'o) which is obtained from the 
mean gel density in a sphere of radius rg around site i. 
The local permeability at site i is then obtained from 
Brinkman's expression. 



hi 



'^.('^o)) 1 



1 - 0i(''o) 



(Al) 



A natural choice for rg is the aerogel correlation length 
^G (or the average cluster size estimated in Ref. ,26. from 



FIG. 16: Histograms of the averaged solid density 1 — <j)i{ro) 
for several values of ro in the 87% (left) and 95% (right) aero- 
gels. 

Finally, we have to introduce the local permeability 
into the evolution equation. This can be done in differ- 
ent ways provided that the conservation of mass is still 
verified and the equilibrium density distribution does not 
change. These conditions are satisfied with the following 
functional form of the evolution equation : 



dt 



J/ 



ij Jij (t) 



(A2) 



where Jij is still defined by Eq. wland the coeffi- 
cients of a permeability matrix which has to be explicited. 
The condition on mass conservation can be written 



as 



4W 



which implies therefore the 

permeability matrix must be symmetric. For simplicity, 
[Ki(fo = Cg) + i^j{rt) = Cg)]/2 



we have chosen Kij 
where ki is defined in Eq. (lAll) 
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It turns out that the fluid density profiles so obtained 
differ very little from those computed with = cf). The 
only deviations occur in the density range correspond- 
ing to the main imbibition front and associated to the 
passage of the liquid in the largest voids of the aerogel. 
These voids fill up more rapidly because the local poros- 



ity in the middle of the cavities reaches high values. As a 
result, the front advances a little bit faster. Overall, the 
differences arc however negligible. One cannot exclude 
that they would increase as time increases, but the pro- 
cedure using local permeabilities is very time-consuming 
and we have not been able to simulate longer times. 
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